#########################################################################
#  File-Name: pobe_appendix_replication.R
#  Date: June 27th 2024
#  Author: Robert Vidigal
#  Purpose: POBE Appendix Replication
#  Out File: pobe_appendix_replication.log
#  Data In: study1dynata.RData; study2fb.RData; study3bovitz.RData
#  Data Out: Analyses, Tables, and Figures for Appendices A-D.
#  Prev file: pobe_replication.R (NECESSARY TO RUN APPENDIX ANALYSES)
#  Status: On-going
#  Machine: Mac Book PRO (Silicon M1 Chip)
#########################################################################

library(lavaan); library(ltm); library(psychometric); 
library(ggplot2); library(ggrepel)

# # -----------------------------------------------------------------------
### APPENDIX A (Gender and DK analyses) 
# # -----------------------------------------------------------------------
source("POBE/scripts/pobe_replication.R")

# # --------------------------------------------------------------------------------------
# Gender differences in giving a “Don’t Know” (Figure 2, Table A.1)
# # --------------------------------------------------------------------------------------

t.test(femalepkMC1DK$spendingdk, malepkMC1DK$spendingdk) # p<.001
t.test(femalepkMC1DK$womendk,    malepkMC1DK$womendk) # p<.001
t.test(femalepkMC1DK$energydk,   malepkMC1DK$energydk) # p=.08
t.test(femalepkMC1DK$merkeldk,   malepkMC1DK$merkeldk) # p<.001
t.test(femalepkMC1DK$vetodk,     malepkMC1DK$vetodk) # p<.001
t.test(femalepkMC1DK$majoritydk, malepkMC1DK$majoritydk) # p<.001

t.test(femalepkTFDK$spendingdk, malepkTFDK$spendingdk) # n.s.
t.test(femalepkTFDK$energydk, malepkTFDK$energydk) # n.s.
t.test(femalepkTFDK$womendk, malepkTFDK$womendk) # n.s.
t.test(femalepkTFDK$vetodk, malepkTFDK$vetodk) # n.s.
t.test(femalepkTFDK$merkeldk, malepkTFDK$merkeldk) # p<.001
t.test(femalepkTFDK$majoritydk, malepkTFDK$majoritydk) # n.s.

# -----------------------------------------------------------------------
### DiD Model 
# -----------------------------------------------------------------------
# Gender-DK -----------------------------------------------------------------
df1$format<-factor(df1$CondM2, levels=c(1,2,5,6), 
                   labels=c("Current", "Current DK", "Certainty", "Certainty DK"))
table(df1$format)

df1$female<-factor(df1$female)
df1$format<-relevel(df1$format, ref = "Current", levels=c(0,1,2,3)); table(df1$format)
did1<-lm(pk~female+format+format*female, data=df1)
summary(did1)

all.effects <- effects::allEffects(mod = did1)
summary(all.effects) 
# Gaps: MC .22; MCDK .87; TF .41; TFDK .45
# DK increases gap among MC, .22 to .87
# DK does not affect the gap as much in TF

# # -----------------------------------------------------------------------
# DIF (STUDY 1)
# # -----------------------------------------------------------------------
### Foreign Aid Spending pk1
### U.S. energy source pk2
### Women SCOTUS pk3
### Angela Merkel office pk4
#### Presidential veto pk5
#### U.S. House majority pk6

# # --------------------------------------------------------------------------------------
### DK0 conditions
# # --------------------------------------------------------------------------------------
# MC
polMC1<-subset(df1, CondM2==1, select=c("spending1", "energy1", "women1", 
                                        "merkel1", "veto1", "majority1", 
                                        "female", "risk")); polMC1<-na.omit(polMC1)
# TF
polTF1<-subset(df1, CondM2==5, select=c("spending3", "energy3", "women3",  
                                        "merkel3", "veto3", "majority3", 
                                        "female", "risk")); polTF1<-na.omit(polTF1)

# Multiple-choice
femaleMC1<-as.character(ifelse(polMC1$female==1,"Female", "Male"))
femaleMC01<-as.numeric(polMC1$female)
polMC1dif<-subset(polMC1, select=c("spending1", "energy1", "women1",  
                                   "merkel1", "veto1", "majority1"))
# True-False with Certainty
femaleTF1<-as.character(ifelse(polTF1$female==1,"Female", "Male"))
femaleTF01<-as.numeric(ifelse(polTF1$female==1,1,0))
polTF1dif<-subset(polTF1, select=c("spending3", "energy3", "women3",  
                                   "merkel3", "veto3", "majority3"))
# Renaming for plotting
polMC1dif <- polMC1dif %>% rename("Foreign Aid Spending"=spending1,
                      "U.S. Source of Energy"=energy1,
                      "Women in SCOTUS"=women1,
                      "Presidential Veto"=veto1,
                      "Angela Merkel's Office"=merkel1,
                      "U.S. House Majority"=majority1)

# Renaming for plotting
polTF1dif <- polTF1dif %>% rename("Foreign Aid Spending"=spending3,
                      "U.S. Source of Energy"=energy3,
                      "Women in SCOTUS"=women3,
                      "Presidential Veto"=veto3,
                      "Angela Merkel's Office"=merkel3,
                      "U.S. House Majority"=majority3)

### CHANGING LATTICE TO PLOT WITH MY DESIRED SETTINGS
# # --------------------------------------------------------------------------------------
mtheme <- lattice::standard.theme("pdf", color=T)
mtheme$strip.border$col ="black"
mtheme$plot.line$lwd <- 3
mtheme$strip.background$col =  "gray90"
mtheme$superpose.line$lwd <- 3
mtheme$superpose.line$col<-c("black", "darkgrey")
mtheme$superpose.symbol$fill<-c("grey30", "lightgrey")

# Multiple Group Model for MC
# # --------------------------------------------------------------------------------------
mGfem1<-multipleGroup(polMC1dif, 1, group = femaleMC1,  
                      invariance=c(colnames(polMC1dif), 'free_means', 'free_var'), SE=T); 
#mGfem1anchor<-multipleGroup(polMC1dif, 1, group = femaleMC1, SE = TRUE,
#                      invariance = c(polMC1dif[,3])); M2(mGfem1anchor)
# Dropdown to allow convergence
dropdownMC <- DIF(mGfem1, c('a1', 'd'), scheme = 'drop')
dropdownMC

# # --------------------------------------------------------------------------------------
# Multiple Group Model for TF
mGfemTF1<-multipleGroup(polTF1dif, 1, group = femaleTF1,  
                        invariance=c(colnames(polTF1dif),'free_means', 'free_var'), SE=T)
# Dropdown to allow convergence
dropdownTF <- DIF(mGfemTF1, c('a1', 'd'), scheme = 'drop')
dropdownTF

# # --------------------------------------------------------------------------------------
# DK1
# # --------------------------------------------------------------------------------------
polMCDK<-subset(df1, CondM2==2, select=c("spending1", "energy1", "women1", 
                                         "veto1", "merkel1", "majority1", 
                                         "female"))
# DK1
# Current
polMCDK1<-subset(polMCDK, select=c("spending1", "energy1", "women1", 
                                   "veto1", "merkel1", "majority1"))
femaleMCDK<-as.character(ifelse(polMCDK$female==1,"Female", "Male"))

# Renaming for plotting
polMCDK1 <- polMCDK1 %>% rename("Foreign Aid Spending"=spending1,
                      "U.S. Source of Energy"=energy1,
                      "Women in SCOTUS"=women1,
                      "Presidential Veto"=veto1,
                      "Angela Merkel's Office"=merkel1,
                      "U.S. House Majority"=majority1)


# Multiple Group Model for MC
# # --------------------------------------------------------------------------------------
mGfem1DK<-multipleGroup(polMCDK1, 1, group = femaleMCDK, SE=T); 
mGfem1DK

# DIF LR TEST, using Wald to ease model convergence
diffemDK1<-DIF(mGfem1DK, which.par = c('a1', 'd'), pairwise = TRUE, Wald=TRUE)
diffemDK1

# # --------------------------------------------------------------------------------------
# # Belief Certainty ---------------------------------------------------------------------
# # --------------------------------------------------------------------------------------
polTFDK<-subset(df1, CondM2==6, select=c("spending3", "energy3", "women3", 
                                         "veto3", "merkel3", "majority3", 
                                         "female"))

polTFDK1<-subset(polTFDK, select=c("spending3", "energy3", "women3", 
                                   "veto3", "merkel3", "majority3"))
femaleTFDK<-as.character(ifelse(polTFDK$female==1,"Female", "Male"))

# Renaming for plotting
polTFDK1 <- polTFDK1 %>% rename("Foreign Aid Spending"=spending3,
                      "U.S. Source of Energy"=energy3,
                      "Women in SCOTUS"=women3,
                      "Presidential Veto"=veto3,
                      "Angela Merkel's Office"=merkel3,
                      "U.S. House Majority"=majority3)

# Multiple Group Model for TF
# # --------------------------------------------------------------------------------------
mGfemTF1DK<-multipleGroup(polTFDK1, 1, group = femaleTFDK, SE=T); 

# DIF LR TEST
# Using Wald to ease model converge
difTFDKfem1<-DIF(mGfemTF1DK, which.par = c('a1', 'd'), pairwise = TRUE, Wald=TRUE)
difTFDKfem1  

# # --------------------------------------------------------------------------------------
### DIF TABLE STUDY 1
# # -----------------------------------------------------------------------
tbldif1<- data.frame("quest"=rep(c("Federal Spending", 
                                   "U.S. Energy Source", 
                                   "Women in SCOTUS", 
                                   "Presidential Veto", 
                                   "Angela Merkel's \n Office", 
                                   "U.S. House \n Majority"), 4),
                     "format"=c(rep('BC',12), rep('MC',12)),
                     "DK"=c(rep("DK0",6), rep("DK1",6), 
                            rep("DK0",6), rep("DK1",6)),
                     "pvalue"=c(dropdownTF$p, difTFDKfem1$p, dropdownMC$p, diffemDK1$p))
tbldif1$format<-factor(tbldif1$format, levels=c("MC", "BC"))
tbldif1$DK<-factor(tbldif1$DK, levels=c("DK1", "DK0"))
tbldif1

# -----------------------------------------------------------------------
# PLOTTING DIF - FIGURE A.1
# -----------------------------------------------------------------------
textdif1<-ggplot(tbldif1, aes(x=format, y=pvalue, color=DK)) + 
  theme_classic() + geom_point() +
  geom_text_repel(label=tbldif1$quest, size = 3, fontface = "bold",
                  min.segment.length = unit(0.05, "npc"), show.legend=F) +
  geom_hline(yintercept=c(.01, .05, .10), linetype="dashed") + 
  scale_y_continuous(limits=c(0,.1), breaks=c(.01, .05, .10)) + 
  ylab("DIF Null Hypothesis Testing \n (p-values)") + 
  xlab("") + theme(legend.position="top", legend.title=element_blank(), 
                   plot.title = element_text(size = 12, face = "bold"), 
                   legend.text=element_text(size=10)) +
  scale_x_discrete(labels=c(paste0("Multiple-choice \n (DK0 N=", nrow(polMC1), 
                                   "; DK1 N=", nrow(polMCDK1), ")"),
                            paste0("Belief Certainty  \n (DK0 N=", nrow(polTF1), 
                                   "; DK1 N=", nrow(polTFDK1), ")"))) + 
  scale_colour_manual(values=c("blue", "black")) + 
  guides(color=guide_legend(ncol=2, override.aes = list(size=4))) +
  theme(axis.text.x = element_text(size=10),  
        axis.title.y = element_text(size=12))

jpeg("POBE/figures/Appendix_Figure_A_1.jpg", 
     units="in", width = 6, height = 5, res=1200)
textdif1
dev.off()

# # -----------------------------------------------------------------------
### DIF NULL HYPOTHESIS TESTING (STUDY 2)
# # -----------------------------------------------------------------------
polMC2<-subset(df2, select=c("know1MC", "know2MC", "know3MC", 
                             "know4MC", "know5MC", "know6MC", "women")); 
polMC2<-na.omit(polMC2)

# Renaming for plotting
# know1 US Term LIMITS
# know2 US Pres VETO
# know3 US House MAJOR
# know4 Women SCOTUS
# know5 US House/Senate Seats up for Election
# know6 UK Prime Minister

# Standard MC
# # --------------------------------------------------------------------------------------
femaleMC2<-as.character(ifelse(polMC2$women==1,"Female", "Male"))
polMC2dif<-subset(polMC2, select=c("know1MC", "know2MC", "know3MC", "know4MC", "know5MC", "know6MC"))
polMC2dif <- polMC2dif %>% rename("U.S. Term Limits"=know1MC,
                      "Presidential Veto"=know2MC,
                      "U.S. House Majority"=know3MC,
                      "Women in SCOTUS"=know4MC,
                      "2022 Election Seats"=know5MC,
                      "U.K. Prime Minister"=know6MC)

# MG
mGfemMC2<-multipleGroup(polMC2dif, 1, group = femaleMC2, SE=T); 

# DIF LR TEST
# Using Wald test and only a1 parameter to allow for model convergence
difMCfem2<-DIF(mGfemMC2, which.par = c('a1'), pairwise = TRUE, Wald=TRUE)
difMCfem2

# Certainty-in-knowledge
# # --------------------------------------------------------------------------------------
polTF2<-subset(df2, select=c("know1TF", "know2TF", "know3TF",  
                             "know4TF", "know5TF", "know6TF", "women")); 
polTF2<-na.omit(polTF2)

# TF
femaleTF2<-as.character(ifelse(polTF2$women==1,"Female", "Male"))
polTF2dif<-subset(polTF2, select=c("know1TF", "know2TF", "know3TF",  
                                   "know4TF", "know5TF", "know6TF"))

polTF2dif <- polTF2dif %>% rename("U.S. Term Limits"=know1TF,
                      "Presidential Veto"=know2TF,
                      "U.S. House Majority"=know3TF,
                      "Women in SCOTUS"=know4TF,
                      "2022 Election Seats"=know5TF,
                      "U.K. Prime Minister"=know6TF)

mGfemTF2<-multipleGroup(polTF2dif, 1, group = femaleTF2, SE=T); 
M2(mGfemTF2, type="C2"); coef(mGfemTF2, IRTpars=F); mGfemTF2 

# DIF LR TEST
# Using Wald test and only a1 parameter to allow for model convergence
difTFfem2<-DIF(mGfemTF2, which.par = c('a1'), pairwise = TRUE, Wald = TRUE); 
difTFfem2

### DIF TABLE STUDY 2
# # -----------------------------------------------------------------------
tbldif2<- data.frame("quest"=rep(c("U.S. Term Limits", 
                                   "Presidential Veto", 
                                   "U.S. House \n Majority", 
                                   "Women in \n SCOTUS", 
                                   "2022 Election Seats", 
                                   "U.K. Prime Minister"), 2),
                     "format"=c(rep('BC',6), rep('MC',6)),
                     "pvalue"=c(difTFfem2$p, difMCfem2$p))
tbldif2$format<-factor(tbldif2$format, levels=c("MC", 'BC'))
tbldif2

# # -----------------------------------------------------------------------
### DIF (STUDY 3)
# # -----------------------------------------------------------------------
### pk1 no term limits
### pk2 prime minister
### pk3 dems house majority
### pk4 3 women supreme court
### pk5 2/3 presidential veto

# MC3
# # --------------------------------------------------------------------------------------
polMC3<-subset(df3, select=c("pk1MC", "pk2MC", "pk3MC", 
                             "pk4MC", "pk5MC", "female")); 
polMC3<-na.omit(polMC3)

polMC3dif<-subset(polMC3, select=c("pk1MC", "pk2MC", "pk3MC", "pk4MC", "pk5MC"))
femaleMC3<-as.character(ifelse(polMC3$female==1,"Female", "Male"))

polMC3dif <- polMC3dif %>% rename("U.S. Term Limits"=pk1MC,
                      "U.K. Prime Minister"=pk2MC,
                      "U.S. House Majority"=pk3MC,
                      "Women in SCOTUS"=pk4MC,
                      "Presidential Veto"=pk5MC)

mGMCfem3<-multipleGroup(polMC3dif, 1, group = femaleMC3, SE=T); M2(mGMCfem3); mGMCfem3 
difMCfem3<-DIF(mGMCfem3, which.par = c('a1', 'd'), pairwise = T, Wald=T); 
difMCfem3

# TF3
# # --------------------------------------------------------------------------------------
polTF3<-subset(df3, select=c("pk1TF", "pk2TF", "pk3TF",  
                             "pk4TF", "pk5TF", "female")); 
polTF3<-na.omit(polTF3)

polTF3dif<-subset(polTF3, select=c("pk1TF", "pk2TF", "pk3TF",  "pk4TF", "pk5TF"))
femaleTF3<-as.character(ifelse(polTF3$female==1,"Female", "Male"))

polTF3dif <- polTF3dif %>% rename("U.S. Term Limits"=pk1TF,
                      "U.K. Prime Minister"=pk2TF,
                      "U.S. House Majority"=pk3TF,
                      "Women in SCOTUS"=pk4TF,
                      "Presidential Veto"=pk5TF)

mGfemTF3<-multipleGroup(polTF3dif, 1, group = femaleTF3, SE=T)
difTFfem3<-DIF(mGfemTF3, which.par = c('a1', 'd'), pairwise = TRUE, Wald=T)
difTFfem3

### DIF TABLE STUDY 3
# # -----------------------------------------------------------------------
tbldif3<- data.frame("quest"=rep(c("U.S. Term Limits", 
                                   "U.K. Prime Minister", 
                                   "U.S. House \n Majority", 
                                   "Women in \n SCOTUS", 
                                   "Presidential Veto"), 2),
                     "format"=c(rep('BC',5), rep('MC',5)),
                     "pvalue"=c(difTFfem3$p, difMCfem3$p))
tbldif3$format<-factor(tbldif3$format, levels=c("MC", 'BC'))
tbldif3

# -----------------------------------------------------------------------
# PLOTTING DIF - FIGURE A.2
# -----------------------------------------------------------------------
tbldif23<-rbind(tbldif2, tbldif3)
tbldif23$study<-c(rep("Study 2",12), rep("Study 3",10))
tbldif23$study<-factor(tbldif23$stud, levels=c("Study 2", 'Study 3'))

textdif23<-ggplot(tbldif23, aes(x=format, y=pvalue, color=study)) + 
  theme_classic() + geom_point() + 
  geom_text_repel(label=tbldif23$quest, size = 3, fontface = "bold",
                  min.segment.length = unit(0.05, "npc")) +
  geom_hline(yintercept=c(.01, .05, .10), linetype="dashed") + 
  scale_y_continuous(limits=c(0,.1), breaks=c(.01, .05, .10)) + 
  ylab("DIF Null Hypothesis Testing \n (p-values)") + 
  xlab("") + theme(legend.position="top", legend.title=element_blank(), 
                   plot.title = element_text(size = 12, face = "bold"), 
                   legend.text=element_text(size=10)) +
  scale_x_discrete(labels=c(paste0("Multiple-choice \n (S2 N=", nrow(polMC2), 
                                   "; S3 N=", nrow(polMC3), ")"), 
                            paste0("Belief Certainty  \n (S2 N=", nrow(polTF2), 
                                   "; S3 N=", nrow(polTF3), ")"))) + 
  scale_colour_manual(values=c("darkorange", "sienna4")) + 
  guides(color=guide_legend(ncol=2, override.aes = list(size=4))) +
  theme(axis.text.x = element_text(size=10),  
        axis.title.y = element_text(size=12))
textdif23

jpeg("POBE/figures/Appendix_Figure_A_2.jpg", 
     units="in", width = 6, height = 5, res=1200)
textdif23
#gridExtra::grid.arrange(textdif1, textdif23, ncol = 1) # together
dev.off()

# # -----------------------------------------------------------------------
### APPENDIX B (Answer look-up analysis) 
# # -----------------------------------------------------------------------

# # -----------------------------------------------------------------------
# LOOK UP REGRESSION MODEL (TABLE B.2)
# # -----------------------------------------------------------------------
# Covariates rescaled 0 to 1
df2$education1<-df2$education/6
df2$income1<-df2$income/12
df2$age1<-df2$age/7
df2$pkTF<-ifelse(df2$pkMC==1,0,1)

summary(lm(know_search ~ pkTF + women + age1 + education + income1 + device_mobile, data = df2))

mobile<-subset(df2, device_mobile==1)
nonmobile<-subset(df2, device_mobile==0)
t.test(mobile$know_search, nonmobile$know_search) # no difference on device type

# -----------------------------------------------------------------------
# -----------------------------------------------------------------------
# -----------------------------------------------------------------------

# # -----------------------------------------------------------------------
### APPENDIX C (VALIDATION AND RELIABILITY) 
# # -----------------------------------------------------------------------

# # --------------------------------------------------------------------------------------
# BINARY MODELS WITHOUT COVARIATES
# # --------------------------------------------------------------------------------------
## TOLERANCE OUTCOME BINARY MODEL
reg_tolMCbinary<-lm(tolerance~pkMC, data=df3); summary(reg_tolMCbinary)
reg_tolTFbinary<-lm(tolerance~pkTF, data=df3); summary(reg_tolTFbinary)

margtolMCbinary<-summary(margins::margins(reg_tolMCbinary, terms=c("pkMC")))
margtolTFbinary<-summary(margins::margins(reg_tolTFbinary, terms=c("pkTF")))

## ENGAGEMENT OUTCOME BINARY MODEL
reg_partMCbinary<-lm(participation~pkMC, data=df3); summary(reg_partMCbinary)
reg_partTFbinary<-lm(participation~pkTF, data=df3); summary(reg_partTFbinary)

margpartMCbinary<-summary(margins::margins(reg_partMCbinary, terms=c("pkMC")))
margpartTFbinary<-summary(margins::margins(reg_partTFbinary, terms=c("pkTF")))

## VACCINATION BELIEFS BINARY MODEL
reg_VACCbeliefsMCbinary<-lm(vaccbeliefs~pkMC, data=df3); summary(reg_VACCbeliefsMCbinary)
reg_VACCbeliefsTFbinary<-lm(vaccbeliefs~pkTF, data=df3); summary(reg_VACCbeliefsTFbinary)

margVACCbeliefsMCbinary<-summary(margins::margins(reg_VACCbeliefsMCbinary, terms=c("pkMC")))
margVACCbeliefsTFbinary<-summary(margins::margins(reg_VACCbeliefsTFbinary, terms=c("pkTF")))

## VACCINE MISINFORMATION BINARY MODEL
reg_VACCmisinfoMCbinary<-lm(vaccmisinfo~pkMC, data=df3); summary(reg_VACCmisinfoMCbinary)
reg_VACCmisinfoTFbinary<-lm(vaccmisinfo~pkTF, data=df3); summary(reg_VACCmisinfoTFbinary)

margVACCmisinfoMCbinary<-summary(margins::margins(reg_VACCmisinfoMCbinary, terms=c("pkMC")))
margVACCmisinfoTFbinary<-summary(margins::margins(reg_VACCmisinfoTFbinary, terms=c("pkTF")))

# # -----------------------------------------------------------------------
### FIGURE - PLOT AVG MARGINAL EFFECTS ACROSS PK TYPES (WITHOUT COVARIATES)
# # -----------------------------------------------------------------------
tblmargBINARY<-data.frame(
  "cond"=c(rep("Political Tolerance",2),rep("Political Engagement",2),
           rep("Vaccine Beliefs",2),rep("Vaccine Misinfo",2)),
                    "knowtype"=c(rep(c("Multiple-choice", "Belief certainty"),4)), 
                    "ame"=c(margtolMCbinary[1,2], margtolTFbinary[1,2],
                            margpartMCbinary[1,2], margpartTFbinary[1,2], 
                            margVACCbeliefsMCbinary[1,2], margVACCbeliefsTFbinary[1,2], 
                            margVACCmisinfoMCbinary[1,2], margVACCmisinfoTFbinary[1,2]), 
                    "SE"=c(margtolMCbinary[1,3], margtolTFbinary[1,3],
                           margpartMCbinary[1,3], margpartTFbinary[1,3], 
                           margVACCbeliefsMCbinary[1,3], margVACCbeliefsTFbinary[1,3], 
                           margVACCmisinfoMCbinary[1,3], margVACCmisinfoTFbinary[1,3])
  )
tblmargBINARY
tblmargBINARY$cond<-factor(tblmarg$cond, levels=c("Political Tolerance", "Political Engagement", "Vaccine Beliefs", "Vaccine Misinfo"), ordered=T)

plotmargBINARY<-ggplot(tblmargBINARY, aes(x=cond, y=ame, shape=knowtype, colour=knowtype)) + 
  geom_errorbar(aes(ymin=ame-1.96*SE, ymax=ame+1.96*SE), width=0, size=1, position=position_dodge(.3), na.rm = T) +
  theme_classic() + geom_point(position=position_dodge(.3), na.rm = T, size=3) +
  theme(legend.position="top", legend.title=element_blank(), 
        plot.title = element_text(size = 12, face = "bold"), 
        legend.text=element_text(size=10)) + 
  geom_hline(yintercept = 0, linetype="dashed", alpha=.3) + 
  theme(axis.text.x = element_text(face='bold', size=10), 
        axis.title.y = element_text(size=12), 
        legend.text = element_text(size=14), axis.text.y = element_text(face='bold', size=12)) + xlab("") + 
  scale_x_discrete(labels=c("Political \n Tolerance", "Political \n Engagement", "Vaccine \n Beliefs", "Vaccine \n Misinformation")) + 
  ylim(-.1,.1) + ylab("Coefficient Estimate") +  guides(col=guide_legend(byrow=T)) + 
  scale_colour_manual(labels = c("Belief certainty", "Multiple-choice"), values=c("black", "grey66"))

jpeg("POBE/figures/Appendix_Figure_C_1.jpg", units="in", width = 6, height = 4, res=1200)
plotmargBINARY # 600x400 #
dev.off()


# # --------------------------------------------------------------------------------------
# FIGURE FOR DISCRETE CHANGE IN MFX
# # --------------------------------------------------------------------------------------

# TOLERANCE # --------------------------------------------------------------------------------------
MCdf3tol<-DAMisc::probci(reg_tolMC, df3, 
               changeX=c("pkMC"), xvals = list(pkMC = seq(0,5,1)),
               returnProbs = T, type = "aveEff")

MCdf3tol<-as.data.frame(MCdf3tol$plot.data)
colnames(MCdf3tol)<-c("pk", "mean", "lower", "upper")

TFdf3tol<-DAMisc::probci(reg_tolTF, df3, 
               changeX=c("pkTF"), xvals = list(pkTF = seq(0,5,.5)),
               returnProbs = T, type = "aveEff")

TFdf3tol<-as.data.frame(TFdf3tol$plot.data)
colnames(TFdf3tol)<-c("pk", "mean", "lower", "upper")

tolplot<-rbind(MCdf3tol, TFdf3tol); tolplot$type<-c(rep(0,6), rep(1,11))
tolplot$type<-c(rep("Multiple-choice", 6), rep("Belief certainty",11))
tolplot$type<-factor(tolplot$type, levels=c("Belief certainty", "Multiple-choice"))
tolplot

#ggeffect(reg_tolMC, terms = c("pkMC"))
#print(ggeffect(reg_tolTF, terms = c("pkTF")), n = Inf) # same result

# PARTICIPATION # --------------------------------------------------------------------------------------
MCdf3part<-DAMisc::probci(reg_partMC, df3, 
               changeX=c("pkMC"), xvals = list(pkMC = seq(0,5,1)),
               returnProbs = T, type = "aveEff")
MCdf3part<-as.data.frame(MCdf3part$plot.data)
colnames(MCdf3part)<-c("pk", "mean", "lower", "upper")

TFdf3part<-DAMisc::probci(reg_partTF, df3, 
               changeX=c("pkTF"), xvals = list(pkTF = seq(0,5,.5)),
               returnProbs = T, type = "aveEff")
TFdf3part<-as.data.frame(TFdf3part$plot.data)
colnames(TFdf3part)<-c("pk", "mean", "lower", "upper")

partplot<-rbind(MCdf3part, TFdf3part);
partplot$type<-c(rep("Multiple-choice", 6), rep("Belief certainty",11))
partplot$type<-factor(partplot$type, levels=c("Belief certainty", "Multiple-choice"))
partplot

# VACCINES BELIEFS # --------------------------------------------------------------------------------------
MCdf3vacc<-DAMisc::probci(reg_VACCbeliefsMC, df3, 
                          changeX=c("pkMC"), xvals = list(pkMC = seq(0,5,1)), 
                          returnProbs = T, type = "aveEff")
MCdf3vacc<-as.data.frame(MCdf3vacc$plot.data)
colnames(MCdf3vacc)<-c("pk", "mean", "lower", "upper")

TFdf3vacc<-DAMisc::probci(reg_VACCbeliefsTF, df3, 
               changeX=c("pkTF"), xvals = list(pkTF = seq(0,5,.5)),
               returnProbs = T, type = "aveEff")
TFdf3vacc<-as.data.frame(TFdf3vacc$plot.data)
colnames(TFdf3vacc)<-c("pk", "mean", "lower", "upper")

vaccplot<-rbind(MCdf3vacc, TFdf3vacc)
vaccplot$type<-c(rep("Multiple-choice", 6), rep("Belief certainty",11))
vaccplot$type<-factor(vaccplot$type, levels=c("Belief certainty", "Multiple-choice"))
vaccplot

# MISINFO # --------------------------------------------------------------------------------------
MCdf3misinfo<-DAMisc::probci(reg_VACCmisinfoMC, df3, 
                  changeX=c("pkMC"), xvals = list(pkMC = seq(0,5,1)), 
                  returnProbs = T, type = "aveEff")
MCdf3misinfo<-as.data.frame(MCdf3misinfo$plot.data)
colnames(MCdf3misinfo)<-c("pk", "mean", "lower", "upper")

TFdf3misinfo<-DAMisc::probci(reg_VACCmisinfoTF, df3, 
               changeX=c("pkTF"), xvals = list(pkTF = seq(0,5,.5)),
               returnProbs = T, type = "aveEff")
TFdf3misinfo<-as.data.frame(TFdf3misinfo$plot.data)
colnames(TFdf3misinfo)<-c("pk", "mean", "lower", "upper")

misinfoplot<-rbind(MCdf3misinfo, TFdf3misinfo)
misinfoplot$type<-c(rep("Multiple-choice", 6), rep("Belief certainty",11))
misinfoplot$type<-factor(misinfoplot$type, levels=c("Belief certainty", "Multiple-choice"))
misinfoplot

# MARGINAL EFFECTS PLOTS # -----------------------------------------------------------------------------
tolplot.out<-ggplot(tolplot, aes(x = pk, y = mean, colour=as.factor(type))) + 
  geom_point(size = 1.5) + 
  geom_errorbar(ymin=tolplot$lower, ymax =tolplot$upper, width = 0.05) + 
  theme_classic() + theme(legend.position="top", legend.title=element_blank()) + 
  xlab("Political Knowledge") + 
  ylab("Marginal Effect on Political Tolerance") + 
  scale_colour_manual(labels = c("Belief certainty", "Multiple-choice"), 
                      values=c("black", "grey66")) +
  scale_y_continuous(limits = c(0,.8), labels = function(x) ifelse(x == 0, "0", x)) + 
  scale_x_continuous(breaks = c(seq(0,5,0.5)), labels = function(x) ifelse(x == 0, "0", x)) + 
  geom_line(aes(group=type))

partplot.out<-ggplot(partplot, aes(x = pk, y = mean, col=as.factor(type))) + 
  geom_point(size = 1.5) + 
  geom_errorbar(ymin=partplot$lower, ymax =partplot$upper, width = 0.05) + 
  theme_classic() + theme(legend.position="top", legend.title=element_blank()) + 
  xlab("Political Knowledge") + 
  ylab("Marginal Effect on Political Participation") + 
  scale_colour_manual(labels = c("Belief certainty", "Multiple-choice"), 
                      values=c("black", "grey66")) +
  scale_y_continuous(limits = c(0,.8), labels = function(x) ifelse(x == 0, "0", x)) + 
  scale_x_continuous(breaks = c(seq(0,5,0.5)), labels = function(x) ifelse(x == 0, "0", x)) + 
  geom_line(aes(group=type))

vaccplot.out<-ggplot(vaccplot, aes(x = pk, y = mean, col=as.factor(type))) + 
  geom_point(size = 1.5) + 
  geom_errorbar(ymin=vaccplot$lower, ymax =vaccplot$upper, width = 0.05) + 
  theme_classic() + theme(legend.position="none", legend.title=element_blank()) + 
  xlab("Political Knowledge") + 
  ylab("Marginal Effect on Vaccine Beliefs") + 
  scale_colour_manual(labels = c("Belief certainty", "Multiple-choice"), 
                      values=c("black", "grey66")) +
  scale_y_continuous(limits = c(0,.8), labels = function(x) ifelse(x == 0, "0", x)) + 
  scale_x_continuous(breaks = c(seq(0,5,0.5)), labels = function(x) ifelse(x == 0, "0", x)) + 
  geom_line(aes(group=type))

misinfoplot.out<-ggplot(misinfoplot, aes(x = pk, y = mean, col=as.factor(type))) + 
  geom_point(size = 1.5) + 
  geom_errorbar(ymin=misinfoplot$lower, ymax =misinfoplot$upper, width = 0.05) + 
  theme_classic() + theme(legend.position="none", legend.title=element_blank()) + 
  xlab("Political Knowledge") + 
  ylab("Marginal Effect on Misinformed beliefs") + 
  scale_colour_manual(labels = c("Belief certainty", "Multiple-choice"), 
                      values=c("black", "grey66")) +
  scale_y_continuous(limits = c(0,.8), labels = function(x) ifelse(x == 0, "0", x)) + 
  scale_x_continuous(breaks = c(seq(0,5,0.5)), labels = function(x) ifelse(x == 0, "0", x)) + 
  geom_line(aes(group=type))

jpeg("POBE/figures/Appendix_Figure_C_2.jpg", 
     units="in", width = 6, height = 6, res=1200)
gridExtra::grid.arrange(tolplot.out, partplot.out, 
                        vaccplot.out, misinfoplot.out, nrow = 2)
dev.off()

# -----------------------------------------------------------------------
# -----------------------------------------------------------------------
# -----------------------------------------------------------------------

# # -----------------------------------------------------------------------
### APPENDIX D
# # -----------------------------------------------------------------------

# TABLES D.1-D.3
round(coef(nomMC1, IRTpars=F,simplify=T)$items[1:6,c(1,5)], 3)
round(coef(nomTF1, IRTpars=F,simplify=T)$items[1:6,], 3)

round(coef(nomMC2, simplify=T)$items[1:6,c(1,5)], 3)
round(coef(nomTF2, simplify=T)$items[1:6,], 3)

round(coef(nomMC3, IRTpars=F, simplify=T)$items[1:5,c(1,5)], 3)
round(coef(nomTF3, IRTpars=F, simplify=T)$items[1:5,], 3)

# # --------------------------------------------------------------------------------------
# OPTION CHARACTERISTIC CURVES
# # --------------------------------------------------------------------------------------

# main script: measurement models section

# -----------------------------------------------------------------------
# -----------------------------------------------------------------------
# -----------------------------------------------------------------------

# # -----------------------------------------------------------------------
### APPENDIX E (Demographics and Knowledge Distributions) 
# # -----------------------------------------------------------------------

# Table 1D - Survey Demographics

# Female
prop.table(table(df1$female)) # 52.9%
prop.table(table(df2$female)) # 50.3%
prop.table(table(df3$women)) # 53.4%

# Democrat
prop.table(table(df1$PID))[2] # 42.5%
prop.table(table(df2$partyid))[3] # 43.2%
prop.table(table(df3$pid))[2] # 47.2% 

# Republican
prop.table(table(df1$PID))[1] # 38.3%
prop.table(table(df2$partyid))[2] # 23.5%
prop.table(table(df3$pid))[1] # 24.6%

# Independent/Other
prop.table(table(df1$PID))[3] + prop.table(table(df1$PID))[4] # 19.2%
prop.table(table(df2$partyid))[1] # 33.3%
prop.table(table(df3$pid))[3] + prop.table(table(df3$pid))[4] # 28.2%

# Mean Age
median(df1$age, na.rm=T) # 1974
median(df2$age, na.rm=T) # 45yo
median(df3$age, na.rm=T) # 3 (35yo)

# Education
median(df1$edu, na.rm=T) # 4
median(df2$edu, na.rm=T) # 3 Some college
median(df3$education, na.rm=T) # 5 College degree

# Income
# df1 Not Asked
median(df2$income, na.rm=T) # 5, 40k
median(df3$income, na.rm=T) # 7, 60k

# # -----------------------------------------------------------------------
### KNOWLEDGE DISTRIBUTION STUDY 1
# # -----------------------------------------------------------------------
require(ggplot2)
pkbar <- ggplot(pkMC1, aes(x=pk)) + geom_bar(color="black", fill="gray66", lwd=0.5) + 
  theme_classic() +
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=14, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 125, by = 25), limits=c(0,125), expand=c(0,0)) +
  scale_x_continuous(name ="(A) Political Knowledge \n Multiple-choice \n M = 3.11 (S.D. 1.44)", expand=c(0.05,0), breaks = seq(from = 0, to = 6, by = 1)) +
  theme(axis.ticks.length = unit(.2,"cm"))

psych::describe(pkMC1$pk) # 3.11; 1.44; median3
var(pkMC1$pk, na.rm=T)

pk3bar <- ggplot(pkTF, aes(x=pk3)) + geom_bar(color="gray66", fill="black", lwd=0.5) + 
  theme_classic() + 
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=14, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 125, by = 25), limits=c(0,125), expand=c(0,0)) +
  scale_x_continuous(name ="(B) Political Knowledge \n Belief Certainty \n M = 2.34 (S.D. 1.12)", 
                     expand=c(0.05,0), limits=c(-0.5,6.5), breaks = seq(from = 0, to = 6, by = .5),
                     labels = function(x) ifelse(x == 0, "0", x)) +
  theme(axis.ticks.length = unit(.2,"cm"))

psych::describe(df1$pk3) # 2.2; 1.24; median2
var(pkTF$pk3, na.rm=T)

t.test(pkMC1$pk, pkTF$pk3) # .76, p<.001

# # --------------------------------------------------------------------------------------
# with DK response option
# # --------------------------------------------------------------------------------------
require(ggplot2)
pkbarDK <- ggplot(pkMC1DK, aes(x=pk)) + geom_bar(color="gray26", fill="gray86", lwd=0.5) + 
  theme_classic() + # geom_bar(data=pkMC1DK, fill = "black", color="black", width=.25, just = -1.5) +
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=14, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 125, by = 25), limits=c(0,125), expand=c(0,0)) +
  scale_x_continuous(name ="(C) Political Knowledge \n Multiple-choice with DK \n M = 2.85 (S.D. 1.49)", expand=c(0.05,0), breaks = seq(from = 0, to = 6, by = 1)) +
  theme(axis.ticks.length = unit(.2,"cm"))


psych::describe(pkMC1DK$pk) # 2.85; 1.49; Median 3
var(pkMC1DK$pk, na.rm=T) # with DK

pk3barDK <- ggplot(pkTFDK, aes(x=pk3)) + geom_bar(color="gray86", fill="gray26", lwd=0.5) + 
  theme_classic() + 
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=14, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 125, by = 25), limits=c(0,125), expand=c(0,0)) +
  scale_x_continuous(name ="(D) Political Knowledge \n Belief Certainty with DK \n M = 2.05 (S.D. 1.34)", 
                     expand=c(0.05,0), limits=c(-0.5,6.5), breaks = seq(from = 0, to = 6, by = .5),
                     labels = function(x) ifelse(x == 0, "0", x)) +
  theme(axis.ticks.length = unit(.2,"cm"))

psych::describe(pkTFDK$pk3) # 2.05; 1.34; Median2
var(pkTFDK$pk3, na.rm=T)

# Figure Study 1 (with and without DK)
jpeg("POBE/figures/Appendix_Figure_E_1.jpg", units="in", width = 12, height = 12, res=1200)
gridExtra::grid.arrange(pkbar, pk3bar, pkbarDK, pk3barDK, ncol = 2) 
dev.off()

t.test(pkMC1DK$pk, pkTFDK$pk3)

# # ----------------------------------------------------------------------
# KNOWLEDGE DISTRIBUTION STUDY 2
# # -----------------------------------------------------------------------
pkMCdf2<-subset(df2, pkMC==1)
pkTFdf2<-subset(df2, pkMC==0)

require(ggplot2)
pkbar2 <- ggplot(pkMCdf2, aes(x=pkMCdf2$knowMCscale)) + geom_bar(color="black", fill="gray", lwd=0.5) + 
  theme_classic() + 
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=10, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 250, by = 50), limits=c(0,250), expand=c(0,0)) +
  scale_x_continuous(name ="(A) Political Knowledge \n Multiple-choice \n M = 2.76 (S.D. 1.55)", expand=c(0.05,0), breaks = seq(from = 0, to = 6, by = 1)) +
  theme(axis.ticks.length = unit(.2,"cm"))

psych::describe(pkMCdf2$knowMCscale)
var(pkMCdf2$knowMCscale, na.rm=T)

pk3bar2 <- ggplot(pkTFdf2, aes(x=knowTFscale)) + geom_bar(color="gray", fill="black", lwd=0.5) + 
  theme_classic() + 
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=10, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 250, by = 50), limits=c(0,250), expand=c(0,0)) +
  scale_x_continuous(name ="(B) Political Knowledge \n with Certainty \n M = 2.83 (S.D. 1.66)", limits=c(-0.5,6.5), breaks = seq(from = 0, to = 6, by = .5),
                     labels = function(x) ifelse(x == 0, "0", x)) +
  theme(axis.ticks.length = unit(.2,"cm"))

psych::describe(pkTFdf2$knowTFscale)
var(pkTFdf2$knowTFscale, na.rm=T)

jpeg("POBE/figures/Appendix_Figure_E_2.jpg", units="in", width = 12, height = 6, res=1200)
gridExtra::grid.arrange(pkbar2, pk3bar2, ncol = 2) # warnings for NAs
dev.off(); rm(pkbar2, pk3bar2) 

# # ----------------------------------------------------------------------
# KNOWLEDGE DISTRIBUTION STUDY 3
# # -----------------------------------------------------------------------
pkMCbar3 <- ggplot(pkMCdf3, aes(x=pkMC)) + geom_bar(color="black", fill="gray", lwd=0.5) + 
  theme_classic() + 
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=10, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 350, by = 50), limits=c(0,350), expand=c(0,0)) +
  scale_x_continuous(name ="(A) Political Knowledge \n Multiple-choice \n M = 3.13 (S.D. 1.21)", expand=c(0.05,0), breaks = seq(from = 0, to = 6, by = 1)) +
  theme(axis.ticks.length = unit(.2,"cm"))

psych::describe(pkMCdf3$pkMC)
var(pkMCdf3$pkMC, na.rm=T)

pkTFbar3 <- ggplot(pkTFdf3, aes(x=pkTF)) + geom_bar(color="gray", fill="black", lwd=0.5) + 
  theme_classic() + 
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) + 
  theme(axis.title.x = element_text(size=18), 
        axis.text.x = element_text(size=10, color="black"),
        axis.title.y = element_text(vjust = 1.5, angle = 90, size = 12), 
        axis.text.y = element_text(size=14, color="black")) +
  scale_y_continuous(name ="Frequency", breaks = seq(from = 0, to = 350, by = 50), limits=c(0,350), expand=c(0,0)) +
  scale_x_continuous(name ="(B) Political Knowledge \n with Certainty \n M = 2.58 (S.D. 1.39)", expand=c(0.05,0), limits=c(-0.5,6), breaks = seq(from = 0, to = 5, by = .5),
                     labels = function(x) ifelse(x == 0, "0", x)) +
  theme(axis.ticks.length = unit(.2,"cm"))

psych::describe(pkTFdf3$pkTF)
var(pkTFdf3$pkTF, na.rm=T)

jpeg("POBE/figures/Appendix_Figure_E_3.jpg", units="in", width = 12, height = 6, res=1200)
gridExtra::grid.arrange(pkMCbar3, pkTFbar3, ncol = 2) # warnings for NAs
dev.off(); rm(pkMCbar3, pkTFbar3) 

# END
print("Code logged successfully")

